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Efficient and robust calculation of femtoscopic correlation functions in spherical 
harmonics directly from the raw pairs measured in heavy-ion collisions 
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We present the formalism for calculating the femtoscopic correlation function directly in spherical 
harmonics. The numerator and denominator are stored as a set of one- dimensional histograms 
representing the spherical harmonic decompositions of each. We present the formalism to calculate 
the correlation function from them directly, without going to any three-dimensional histogram. We 
discuss the practical implementation of the method and we provide an example of its use. We also 
discuss the stability of the method in the presence of angular holes in the underlying data (e.g. from 
experimental acceptance). 

PACS numbers: PACS numbers: 25.75.-q, 25.75.Gz 



Spherical harmonics are one of the most commonly 
used mathematical tools for the analysis of experimental 
data. For example, geopotcntial models of the Earth's 
gravitational field are matched to experimental data up 
to harmonic order £ rnax = 70 Measurements of the 
cosmic microwave background is expanded to by the Cos- 
mic Background Explorer (COBE) and successor ex- 
periments extend further to an impressive £ max ~ 1500 
|3j . In the femtoscopic measurements of particle emitting 
sources in heavy-ion collisions, one also finds it useful to 
expand the measured correlations and extracted sources 
in spherical harmonics 0, 0] ■ 

In nearly all applications of spherical harmonics to 
data analysis and reduction, one is faced with theprob- 
lem of "holes in the data" - i.e. sampling bias [6j. In 
the COBE analysis mentioned above, this bias occurs 
because the Milky Way masks a sizable solid angle of the 
sky. When constructing potential maps of the earth, the 
problem is even more severe as one uses strips of data 
obtained from various satellite and balloon-borne exper- 
iments to derive the map Q- In heavy-ion collisions, the 
sampling bias most often arises because the detector ac- 
ceptance does not span all of phase-space. 

Unlike other applications, the femtoscopic correlation 
functions are actually ratios of two single particle distri- 
butions: the true pair distribution (the numerator) and 
the mixed pair distribution (the denominator). Thus, ex- 
panding the correlation function in spherical harmonics 
is a more involved than expanding the underlying single 
particle distributions. It is often impractical to simply 
bin the two singles distributions in 3D histograms and 
make a ratio since one does not have a meaningful num- 
ber of pairs in each bin, either due to the statistics of 
particle production or due to detector acceptances. Fur- 
thermore, since the end result is a spherical harmonic de- 
composition of the correlation, one must have absurdly 
high harmonics in one's decomposition to resolve high- 



momentum bins which have small angular extent and 
marginal statistics. 

Rather than pursue this, we adopt an alternate ap- 
proach: we construct the raw pair distributions directly 
in spherical harmonics, then we extract the correlation 
function, also in spherical harmonics, by viewing the ra- 
tio as an inverse problem. A feature of this approach is 
that one preserves the full cross-£, m data covariance that 
is currently ignored when imaging 3D correlations 

We now outline this paper. In the first section, we de- 
tail how we expand the numerator and denominator dis- 
tributions in spherical harmonics and how we pack them 
into the vectors and matrices for further manipulation. 
In the second section, we describe how to compute the 
correlation in spherical harmonics using these harmonic 
expansions. In the last section, we demonstrate the tech- 
nique in realistic examples. Further, we show how this 
technique is insensitive to the sampling bias imposed by 
sizable holes in pair acceptance. 



CONSTRUCTING THE NUMERATOR AND 
DENOMINATOR 

The femtoscopic correlation function, C(q), is defined 
as a ratio of the probability to observe a correlated pair of 
particles at a given relative momentum in the same event 
(numerator), T(q), to the probability to observe such a 
pair in an uncorrelated state (denominator), M(q): 



T(q) = C(q)M(q). 



(1) 



The uncorrelated distribution is usually obtained by mix- 
ing particles from different events. Here the relative mo- 
mentum is given in the pair center of mass frame as 
q = |(Pi — P2) = k* for particles 1 and 2. Here k* is 
the notation commonly used for femtoscopic correlations 
for pairs of two different types of particles, and will be 
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Condition 



Relation 



1 Distribution is real Tt m = T/_ m 

2 x -> -x symmetry T im (x,y,z) = {-l) m T^ m (-x,y, z) 

3 y -» symmetry T^fz, y, z) = T/ m (a;, -y, z) 

4 z -> -2 symmetry T lm (x, y, z) = (~l) l+m T lm {x, y, -z) 

5 r — » -r symmetry T im (x,y, z) = (~l) e T em (-x, -y, -z) 

TABLE I: Conditions on the distribution imply relations be- 
tween the different terms in the Yi m expansion. Condition 
1 is always valid for the true and mixed pair distributions 
and the correlation function. Condition 5 is always valid for 
like pair correlation and the corresponding pairs distributions. 
By exploiting symmetries in 2-5, we can reduce the number 
of components in our data vectors. 



(if qi in radial bin n only) : 
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(5) 



Taking the diagonal elements (i.e. im = £'m'), we find 
(uncorrelated) uncertainties of 



AT J f m (q n ) sa A 2 Te m e m (q n ). 



(6) 



used interchangeably with q (traditionally used in femto- 
scopic correlations of pairs of identical partilccs) later in 
the paper. We expand the numerator, denominator and 
correlation function in spherical harmonics, e.g.: 



T(q) 



7T J~] Tj m (q) Yj m (f2 q ) . 

tm 



(2) 



We can compute the pairs distributions directly in 
spherical harmonics by observing that 



Tt m {q) 



1 



dn 4 r(q)y; m (n 4 ). 



(3) 



can be built up by summing over the pairs, which is es- 
sentially a Monte-Carlo integration process: 
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otherwise. 



(4) 



Here, N is the number of pairs in the spectrum. In our 
approach a pair is added to Ti m (q n ) (Mt m (q n )) in the 
following way. First we calculate the relative momentum 
q and decompose it into |q| (which determines the ID 
g-bin number) , # qi and <^ qi . Having the angles, we can 
calculate the spherical harmonics functions, usually up 
to some limiting value of £. Then the pair is added to all 
histograms of the corresponding function, with weights 
equal to the respective Yg^s. We could eliminate many 
of the components of Te m (q n ) and Mg m {q n ) by taking 
advantage of the symmetries in Table I, but we have 
chosen to keep all components so that we can perform 
cross-checks of our work. This does introduce complica- 
tions when performing some matrix manipulations, as we 
discuss in the following section. 

The covariance can be built in a similiar way by noting 



In contrast with our approach, in the traditional rep- 
resentation both the numerator and denominator were 
stored as 3D histograms, using either the Bertsch-Pratt 
coordinates [f| in Cartesian form q ou t, qside an d qiong, or 
in spherical: |q|, cos(# q ), </> q . In the traditional represen- 
tation the numerator (denominator) is a 3D histogram, 
and each signal (background) pair is added to exactly 
one bin of this histogram with weight 1.0. This represen- 
tation has several disadvantages: one needs significant 
statistics to have a meaningful number of pairs in each 
bin, single-particle momentum acceptance can result in 
"holes" or empty bins in two-particle, or relative mo- 
mentum space, and the binning corrections need to be 
applied. Also going to higher moments in the decompo- 
sition requires larger number of bins on the </> q and 6* q 
direction. 

In practice, we store both the real and imaginary parts 
of the numerator and the denominator as an array of one- 
dimensional histograms in q = |q|, as seen in Eq. |7|. 
which is represented as the vector T ? : 



Too(q) 
T W (q) 

RTii(?) 

9Tii(g) 

T 20 (q) 
^T 21 {q) 
%T 21 (q) 
®T 22 {q) 
%T 22 (q) 

T 3 o{q) 



(7) 



There is also a corresponding covariance matrix A 2 T 
which is a two-dimensional matrix for each g-bin and is 
packed in an analogous fashion. 



?! 
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CALCULATING THE CORRELATION 
FUNCTION 

Since we have expanded the pair distributions and the 
correlation in spherical harmonics, we have 



In order to calculate the correlation function, we view 
the problem as an inverse problem. To solve it, one needs 
to minimize the x ■ 

(T, - M q ■ C q ) T ■ (A 2 ^)- 1 • (T, - M, • C,). (10) 

The formula uses the full covariance matrix in the true 
distribution, but not in the mixed pair distribution. Be- 
cause the mixed pair distribution is constructed by pairs 
from different events, it is not limited by statistics and 
can be computed to arbitrarily high precision, making 
the uncertainties negligible for our purposes. The prob- 
lem of minimizing the \ 2 is identical to the one posed by 
the imaging procedure in Ref. Q and the solution is well 
known: 

C q = A 2 C q ■ M.g ■ (A 2 T g ) _1 • T q . (11) 

where the covariance is also calculated: 

A 2 C q = (M£ • (A 2 ^)- 1 ■ M,)- 1 . (12) 

The uncorrelated uncertainties of the correlation are just 
the square root of the trace of the covariance A 2 C g : 
AC g = ^TrA 2 C q . Written out, we are simply tak- 
ing the diagonal elements (i.e. tm = £'m') and finding 
(uncorrelated) uncertainties of 

AQm(?n) ~ \/A 2 CW m (g„). (13) 

With the error propagation done in this way, one gets the 
cross-^m correlations in the covariance matrix by virtue 
of the cross-£m correlations built into the M g and A 2 T g 
matrices. 

We compute the correlation function according to 
Eq. (jlip . As this involves several matrix inversions, it 
is important to make sure that the matrix determinant 
is not zero. That means that one cannot use £m combina- 
tions that arc a linear combination of other Em's. As all 



With the packing in Eq. ((7|) , equation JH]) can be written 
very compactly: T q — M 9 ■ C q . Eq. ([8]) gives us a way 
to compute C'e m (q) directly from the pair distributions 
expanded in spherical harmonics. 



(9) 

I 

of the pair distributions are real, we cannot keep m > 
and m < components at the same time (see Table [J) . 
Therefore we adopt the convention to only use positive 
components (including m = component). Therefore we 
remove all functions with negative m from both T q and 
M 9 when solving Eq. (fTTjl . One can add the missing neg- 
ative m components to the correlation function by multi- 
plying the positive m values by the appropriate factor of 
(— l) l+m . For consistency one may repeat the procedure, 
this time removing the positive m components from T q 
and M g and the obtained results should be identical. 

We remind the reader that each q bin is independent 
in T g , M g and C q . Therefore solving Eq. (fTTj) can be 
done for each q bin independently. The starting points 
then are not vectors of functions, but simply vectors of 
real numbers T q and M g and the result is also a vector 
of real numbers C q . A covariance matrix (A 2 T) g is in 
this case a 2D matrix of real numbers, as is the resulting 
correlation function covariance matrix (A 2 C) 9 . Solution 
of Eq. (TTT]) then reduces to a problem of solving a set 
of linear equations, for which many standard numerical 
algorithms exist. The procedure is repeated for each q 
bin and the C q vector is filled in steps. 



REALISTIC EXAMPLES 

We have tested and applied the formalism to the con- 
struction of two-particle correlation functions for iden- 
tical and non-identical particles in relativistic heavy-ion 
collisions. Our tests include a study of the robustness of 
our approach in the presence of 9 — 4> acceptance holes 
in relative momentum, a common occurance in exper- 
iments including the STAR experiment of which one of 
the authors (Kisicl) is a collaboration member. Below we 
present some tests of the method with a realistic model. 



Ti m (q) = 2_j M t'm'{q)Ci>" m "{q) 
t'm' l"m" 

x / ^(fi^^y^^aa) (8) 

Jin 

= ^■■if"m"(g)^<"m"(g) - 1 

Here the Mj'matrix is written in terms of Wigner 3-j symbols as: 

M lml „ m „ = £ Me m ,(q)(-l) m VW + W' + W' + 1) ( t l ' e ")( 1 e , ) ■ 
JT^i 1 M — m m m I 
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FIG. 1: (Color online) Numerator of an example tv + K + correlation function binned directly in spherical harmonics, as a function 
of the first particle's momentum in the pair rest frame k*. Panel a) shows 1 = 0, panels b) and c) show I = 1 components, 
center panels d), e) an f) shows 1 = 2 components, right panels g), h), i) and j) show £ = 3 components. Red open circles 
represent real part of the decomposition, blue closed triangles show imaginary part. 



Example correlation functions 

Two particles are femtoscopically correlated if they 
have a small relative momentum k* in the pairs' rest 
frame (for identical particles we use q = k*). If they are 
not identical and have different masses, they must have 
different momenta in the laboratory frame. In such case 
it can happen that due to specific momentum acceptance 
of the experiment, pairs with specific values of k* = |k*| 
and certain combinations of polar and azimuthal compo- 
nents of k* cannot be measured. In terms of the spherical 
harmonic representation, this results in a hole in the pair 
acceptance for certain regions in k*, (f> and cos 8. This 
is observed e.g. for pion-kaon pairs in the STAR experi- 
ment. Such a hole presents a methodological problem for 
traditional methods of decomposing the correlation func- 
tion in spherical harmonics as they rely on the existence 
of certain symmetries in pair distributions. In particular, 
they assume that the multiplicity of pairs with a given 
k* ut is equal to the multiplicity of pairs with —k* ut . For 
non-identical particles there is no such symmetry While 
it is certainly possible to improve the existing methods 



and to remove this dependence, we propose to move to 
the more advanced decomposition method presented in 
this paper and bypass the problem altogether. In our 
method, this hole is reflected in both the numerator and 
denominator by a lower number of pairs contributing at 
some bins of k*. 

Examples of the numerator of the correlation func- 
tion, binned directly in spherical harmonics, are shown 
in Fig. [TJ In this case, the distributions result from a 
simulation of the tt + K + correlation function in the Ther- 
minator model using the STAR detector acceptance [§]. 
The acceptance is symmetric with respect to cos 6 so 
the 1, m = 1,0 component vanishes [f|. All the imag- 
inary components vanish. Also the i, m = 2, 1 as well 
as £, m = 3, and £, m — 3, 2 vanish due to polar angle 
symmetry. Apart from that, the numerator shows non- 
trivial structure both as a function of k* and 4>. The 
synergy between spherical harmonic decomposition and 
femtoscopic correlation function is nicely illustrated in 
this plot. A full 3D information, which in traditional 3D 
implementation would require tens of thousands of bins 
to store, is reduced to a few ID histograms. Out of these 
only a select few carry important information, while oth- 
ers conveniently vanish due to the intrinsic symmetries 
of the pair distribution. The significance of the compo- 
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nents diminishes with growing £, ensuring that cutting 
the decomposition at some £ m ax should not distort the 
function. 

Having in mind that the underlying numerator has a 
non-trivial structure both in k* and <j> it is interesting 
to see how the correlation function itself, calculated with 
the method above, behaves. It is shown in Fig. [2j Again, 
the imaginary components all vanish, as they should. 
The Coo component shows the expected behavior com- 
ing from a Coulomb repulsion of same charge pion and 
kaon. The C20 and C22 components show small devia- 
tions from zero, which signals the fact that the size of 
the underlying system is not the same in the out, side 
and long directions. Also the Cxi component deviates 
from zero significantly, a signature of the average emis- 
sion point asymmetry between pions and kaons. In sum- 
mary the example confirms several important points: (a) 
the correlation function can be calculated via the direct 
Yi m method, (b) The important physics signals in Coo, 
Cli, C20 and C22, are preserved, (c) Other components 
of the correlation function vanish, as they should - an 
important cross-check of the method. 



Limits of applicability in presence of an acceptance 

hole 



Our results are in contrast to what would happen in 
the traditional approach of expanding the correlation in 
spherical harmonics after making the ratio of 3D his- 
tograms. If there is poor statistics due to a gap in accep- 
tance, then one will need a large number of spherical mo- 
ments to capture the purely statistical fluctations present 
in the poorly populated high-q bins. What is more insid- 
ious, because the correlation is a ratio, the structure in 
poorly determined 3D bins appear to "cancel out" even 
when the poorly resolved data should not cancel out. 
Rather, the poor statistics should give rise to large un- 
certainties and not contribute to the spherical harmonic 
expansion (which is what happens in our method). 

In Fig. |H we show the values of the analytical fit to the 
spherical harmonics vs. the hole size. The lines are the 
"input radii" from the Thcrminator model. As one sees, 
the fit results are very stable and moreover in reasonable 
agreement with the input values. Note: exact agreement 
between the input and extracted radii can not be ex- 
pected because the collective motion in the Thcrminator 
model shrinks the effective homogeneity length seen by 
the pairs and hence the correlation radii. Even the very 
large hole of 37r/2 (only a quarter of acceptance remain- 
ing!) our method seems to preserve all the relevant com- 
ponents, provided that enough statistics remains outside 
the hole region. 



To attempt to determine the practical limits of the 
technique, we have performed a test. First, we calcu- 
late the correlation function for identical pions using the 
Thcrminator model. We chose to use identical neutral 
pion pairs for the calculation simply because final state 
interactions do not distort the correlation appreciably, 
meaning that the correlation shape can be characterized 
simply by the correlation radii R S ide, Rout and Ri ong - 
The source size has been set to reasonable values in the 
longitudinally co- moving system (of ~ 3 — 4 fm). The 
first calculation does not have any acceptance holes. We 
then repeat the calculation, introducing an artificial hole 
in the acceptance by removing both from the numerator 
and denominator all pairs within the hole. The hole is 
at midrapidity (small cos (9), small q (from 0.01 to 0.05 
GeV/c) and with varying width in <j> (from to 3ir/2). 
The results are shown in Fig. [3] 

One can see that introducing the hole had no influence 
on the extracted correlation function within statistical 
errors. Indeed, the dominant effect of the acceptance 
hole has been to decrease statistics, increasing statistical 
scatter and the corresponding uncertainty. To further 
make this point, we show the analytical prediction for 
how the spherical harmonics should look like for these 
sizes in the black dashed lines. As one can see, all points 
follow the lines perfectly. 



Experimental corrections 

In order to be useful, the procedure for calculating 
the correlation function directly in spherical harmonics 
should allow for the application of the standard experi- 
mental corrections. Here we briefly describe how this can 
be done. 

Experimental resolution for two-particle reconstruc- 
tion and identification is usually dominated by two issues: 
track merging (where two tracks in the detector are re- 
constructed as one) and track splitting (where a single 
track is mistakenly reconstructed as two). These have 
non-trivial dependence on both the particle momenta as 
well as their trajectory in the detector. This is usually 
corrected for by assigning a weight to each pair, based 
on the detailed detector simulation. Such weighting can 
be incorporated in the procedure in a straightforward 
way. When filling the numerator and the denominator 
with pairs, one simply fills it with the appropriate weight. 
Mathematically it amounts to modifying Eq. (j4]) by mul- 
tiplying the i^C^gJ by an additional weight W, coming 
from the above mentioned correction. 

Another common issue is the particle purity, namely 
the fraction of pairs in the sample that should be treated 
as correlated. A pair may be not correlated if one of the 
particles is misidentified or if at least one of the particles 
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FIG. 2: (Color online) The example correlation function binned directly in spherical harmonics. The red open circles are the 
real parts of each term and the filled blue squares are the imaginary part of each term. The spherical moments of the numerator 
spectra are shown on panels as in Fig. [T] 



comes from a weak decay. The experiment should be able 
to estimate the the average purity of pairs P, which can 
be (and usually is) a function of particles' momenta, and 
therefore also of the pair relative momentum q. We use 
the traditional formula: 

Ccorr(cd= Cme p { ^~ 1 +l. (W) 

From the correlation function C, we can obtain the cor- 
relation effect R = C — 1. In spherical harmonic rep- 
resentation this only modifies the £ = 0, m = com- 
ponent: i?oo = Coo ~ 1 5 while others remain the same: 
Rim = Gi m . Then Eq. (fTi|) simplifies to: 

D I „\ -^meas(q) k n 
Rcorr(<l) = p , s ■ (15) 

We immediately note that it is equivalent to calculating 
the correlation function from the numerator and denomi- 
nator. Therefore it is enough to express purity P directly 
in spherical harmonics and treat it as denominator, take 
the measured correlation function and treat it as numer- 
ator and finally apply the mathematical formalism de- 
scribed in this work to obtain the correlation function 
corrected for purity. 



APPLICABILITY 

The method presented in this paper has been successful 
applied to the femtoscopic correlation functions in heavy- 
ion collisions. It should be possible to apply it to other 
fields as well, however one has to take into account limits 
of the method applicability. 

The basic formula ([5]) is strictly correct mathematically 
only if one uses an infinite number of £, m components for 
all the functions T q , M g and C q . In practical application 
one needs to limit oneself to a specific value of I. This 
is only allowed if the higher ^-moments are negligible. 
The femtoscopic correlation function is very well suited 
to the method because the intrinsic symmetries of the 
pair distributions limit the relevant I components to a 
practical maximum of 6. Most important information is 
contained in I = 0, £ = 1 and £ = 2 components. 

We have shown that the method remains stable for any 
reasonable acceptance hole in <f> region. It is also clear 
that the method will start breaking down only for re- 
ally small values of (f> — 9 acceptance and in the extreme 
case of the "hole" taking up the whole <f> — 9 acceptance 
the method will simply break down due to lack of data. 
The method in effect interpolates the correlation func- 
tion in the region where there is no data by assuming 
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FIG. 3: (Color online) The simulation of the acceptance hole, made at midrapidity (small cos 9) at small q = 2k* (from 0.01 
to 0.05 GeV/c) and with varying width in <j> (from to 3n/2). The plot shows the correlation in spherical harmonics for: ideal 
case with no hole (black diamonds), "small hole" n/6 (red circles), "sizeable hole" ir/2 (green squares), and "huge hole" 3tt/2 
(blue stars). The black dashed lines show the analytical prediction for how the spherical harmonics should look like for these 
sizes. 
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FIG. 4: (Color online) The values of the analytical fit to the 
spherical harmonics vs. the hole size. The solid lines are 
the "input radii" from the Therminator calculations and the 
dashed lines/symbols are our fits. 



certain symmetries in the underlying pair distribution. 
Using the method described here for femtoscopy, accep- 
tance holes are "irrelevant" - any reasonable femtoscopic 
measurement will have a large enough acceptance to be 
insensitive to the holes. 
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